outliers,
предположили, что выбросы это значения, которые отклоняются более чем на
1.5 межквартильного размаха от первого и третьего квартилей.Новые библиотеки:
plotly: визуализация данных (построение
графиков)
patchwork: объединение графиков и более гибкая
настройки визуализации результатов ggplot2
effects, sjPlot: проверка условий
применимости смешанных линейных моделей
purrr: упрощает сложные манипулцяии с данными, с
помощью функции map в пакете возможно значительно сократить
описание функции и переменных
Данная работа посвящена анализу сортов льна, выращиваемым в течение 5 лет на 2 разных локациях: на Кубани и в Липецке. Начнём анализ с изучения переменных, по которым была собрана информация:
## 'data.frame': 2943 obs. of 16 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ leaf_shape : chr "round" "round" "round" "round" ...
## $ maturation_group : int 5 6 5 3 5 3 5 5 6 5 ...
## $ lodging_type : chr "yes" "yes" "no" "no" ...
## $ growth_type : chr "indeterminant" "indeterminant" "indeterminant" "semi_determinant" ...
## $ flowering_group : num 5 5 5 4 5 2 4 5 5 5 ...
## $ pubescence_colour: chr "light_gray" "light_gray" "light_gray" "light_gray" ...
## $ corolla_colour : chr "white" "purple" "purple" "purple" ...
## $ origin : chr "Uzbekistan" "India" "China" "China" ...
## $ productivity : num NA NA NA NA NA NA NA NA NA NA ...
## $ vegetation_period: int NA NA NA NA NA NA NA NA NA NA ...
## $ protein_content : num NA NA NA NA NA NA NA NA NA NA ...
## $ oil_content : num NA NA NA NA NA NA NA NA NA NA ...
## $ site : chr "kub" "kub" "kub" "kub" ...
## $ year : int 2017 2017 2017 2017 2017 2017 2017 2017 2017 2017 ...
Заметим, что в исходном датасете есть пропущенные значения и категориальные переменные, которые для дальнейшей работы необходимо будет перевести в факторные.
| Переменная | Описание | Вариации |
|---|---|---|
| id | Номер сорта | 300 сортов |
| leaf_shape | Форма листа | lanceolate/round |
| maturation_group | Группа созревания | от 1 (раннеспелый) до 6 (позднеспелый) |
| lodging_type | Тип полегания | no/leaning/yes |
| growth_type | Тип развития | determinant/indeterminant/semi_determinant |
| flowering_group | Группа цветения | от 1.00 (рано) до 5 (поздно), шаг 0.5 |
| pubescence_colour | Цвет опушения | gray/light_gray/light_tawny/tawny |
| corolla_colour | Цвет венчика | purple/white |
| origin | Происхождение | Страна |
| productivity | Продуктивность | г/м2 |
| vegetation_period | Период от всхода до сбора урожая | дни |
| protein_content | Содержание белка в семенах | 1-100% |
| oil_content | Содержание масла в семенах | 1-100% |
| site | Место выращивания | kub (Кубань), lip (Липецк) |
| year | Год сбора | 2017-2021 |
В выборке встречаются сорта родом из 31 страны.
## [1] "Austria" "Belarus" "China" "Czech" "Estonia"
## [6] "France" "Georgia" "Germany" "Hungary" "India"
## [11] "Italy" "Japan" "Kanaues" "Kazakhstan" "Latvian"
## [16] "Lithuania" "Moldova" "Netherlands" "NorthKorea" "Poland"
## [21] "Portugal" "Romania" "Russia" "Serbia" "Slovakia"
## [26] "South_Korea" "Sweden" "Ukraine" "USA" "Uzbekistan"
## [31] "Yugoslavia"
Переменные leaf_shape, lodging_type,
growth_type, pubescence_colour,
corolla_colour, origin, site
переведём в факторные.
Maturation_group и flowering_group — в
упорядоченные факторные.
Productivity, vegetation_period,
protein_content, oil_content - оставим в
первозданном виде.
data <- data %>%
mutate(across(c(
leaf_shape, maturation_group, lodging_type, growth_type, flowering_group,
pubescence_colour, corolla_colour, site, year, origin), as.factor)) %>%
mutate(
maturation_group = factor(maturation_group,
levels = c(1, 2, 3, 4, 5, 6),
ordered = TRUE),
flowering_group = factor(flowering_group,
levels = c(1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5),
ordered = TRUE)
)Рассмотрим данные с помощью функции summary
## X id leaf_shape maturation_group
## Min. : 1.0 Min. : 1.0 lanceolate: 225 1: 117
## 1st Qu.: 736.5 1st Qu.: 84.0 round :2718 2: 468
## Median :1472.0 Median :167.0 3:1017
## Mean :1472.0 Mean :166.4 4: 828
## 3rd Qu.:2207.5 3rd Qu.:249.0 5: 378
## Max. :2943.0 Max. :330.0 6: 135
##
## lodging_type growth_type flowering_group pubescence_colour
## leaning: 162 determinant : 936 3 :1449 gray : 90
## no :2655 indeterminant : 729 5 : 414 light_gray :1701
## yes : 126 semi_determinant:1278 4 : 387 light_tawny:1143
## 3.5 : 234 tawny : 9
## 4.5 : 198
## 2 : 162
## (Other): 99
## corolla_colour origin productivity vegetation_period
## purple:2331 Russia :1557 Min. : 9.9 Min. : 72.0
## white : 612 Ukraine: 189 1st Qu.:128.0 1st Qu.: 96.0
## France : 153 Median :173.0 Median :104.0
## Kanaues: 126 Mean :170.2 Mean :103.9
## China : 117 3rd Qu.:210.0 3rd Qu.:113.0
## (Other): 693 Max. :362.0 Max. :137.0
## NA's : 108 NA's :1543 NA's :1533
## protein_content oil_content site year
## Min. :33.00 Min. :14.30 kub:1635 2017:327
## 1st Qu.:39.50 1st Qu.:20.00 lip:1308 2018:654
## Median :41.20 Median :21.20 2019:654
## Mean :41.07 Mean :21.05 2020:654
## 3rd Qu.:42.60 3rd Qu.:22.30 2021:654
## Max. :51.50 Max. :25.10
## NA's :1686 NA's :1686
Поскольку в датасете довольно много категориальных описательных переменных мы проверили их на наличие опечаток и ничего не обнаружили.
## [1] lanceolate round
## Levels: lanceolate round
## [1] determinant indeterminant semi_determinant
## Levels: determinant indeterminant semi_determinant
## [1] leaning no yes
## Levels: leaning no yes
## [1] gray light_gray light_tawny tawny
## Levels: gray light_gray light_tawny tawny
## [1] purple white
## Levels: purple white
## [1] kub lip
## Levels: kub lip
Всего в датасете у нас 6556 пропущенных значений. На первый взгляд их довольно много. Подробнее изучим в каких столбцах они встречаются:
origin содержит 108 пропущенных значений.productivity содержит 1543 пропущенных значений.vegetation_period содержит 1533 пропущенных
значений.protein_content содержит 1686 пропущенных
значений.oil_content содержит 1686 пропущенных значений.При удалении всех пропущенных значений объем данных для дальнейшего анализа уменьшается в два раза, но при этом еще и полностью пропадают данные о познеспелых сортах (6), поскольку по ним нет измерений по количественным признакам. Кроме того, мы заметили, что за 2017 год в Липецке не проводилось измерений:
##
## kub lip
## 2017 327 0
## 2018 327 327
## 2019 327 327
## 2020 327 327
## 2021 327 327
Мы можем позволить себе удалить все пропущенные значения, поскольку данных будет всё ещё достаточно для проведения статистического анализа:
data_without_na <- data %>%
filter(if_all(everything(), ~ !is.na(.) & . != "" & trimws(.) != ""))
summary(data_without_na)## X id leaf_shape maturation_group
## Min. : 141.0 Min. : 1.0 lanceolate: 105 1: 54
## 1st Qu.: 467.5 1st Qu.:117.0 round :1090 2:260
## Median : 793.0 Median :179.0 3:510
## Mean :1078.5 Mean :177.0 4:259
## 3rd Qu.:1705.5 3rd Qu.:239.5 5:112
## Max. :2826.0 Max. :330.0 6: 0
##
## lodging_type growth_type flowering_group pubescence_colour
## leaning: 69 determinant :346 3 :693 gray : 38
## no :1098 indeterminant :282 4 :129 light_gray :685
## yes : 28 semi_determinant:567 2 : 99 light_tawny:470
## 3.5 : 91 tawny : 2
## 5 : 76
## 4.5 : 59
## (Other): 48
## corolla_colour origin productivity vegetation_period protein_content
## purple:942 Russia :687 Min. : 12.0 Min. : 72 Min. :33.00
## white :253 France : 85 1st Qu.:132.0 1st Qu.: 96 1st Qu.:39.50
## Ukraine: 81 Median :177.0 Median :104 Median :41.20
## Kanaues: 80 Mean :173.7 Mean :104 Mean :41.02
## Austria: 45 3rd Qu.:213.0 3rd Qu.:113 3rd Qu.:42.60
## Serbia : 36 Max. :362.0 Max. :137 Max. :51.50
## (Other):181
## oil_content site year
## Min. :14.30 kub:844 2017:173
## 1st Qu.:20.00 lip:351 2018:401
## Median :21.20 2019:351
## Mean :21.06 2020:161
## 3rd Qu.:22.40 2021:109
## Max. :25.10
##
Проанализируем выбросы в численных переменных
productivity, protein_content,
oil_content. Стоит учесть, что это непростая задача,
поскольку в выброке 300 разных сортов, обладающих уникальными
признаками, потому выборка гетерогенна и сорта высоко/низкопродуктивные
могут быть приняты как выброс.Поэтому данные для удаления стоит
анализировать в совокупности с остальными признаками.
Будем искать с помощью пакета outliers. Предположим, что
выбросы это значения, которые отклоняются более чем на 1.5
межквартильного размаха от первого и третьего квартилей.
outliers <- boxplot.stats(data_without_na$productivity)$out
ggplot(data_without_na, aes(x = maturation_group, y = productivity, group = maturation_group)) +
geom_boxplot(
outlier.shape = NA,
color = "black",
lwd = 1.0
) +
geom_jitter(
aes(color = ifelse(productivity %in% outliers, "outlier", "normal")),
width = 0.2,
height = 0,
size = 1
) +
scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +
theme_minimal(base_size = 14) +
labs(
x = "Группа созревания",
y = "Продуктивность",
color = "Тип данных"
)
2. Содержание белка и содержание масла
Можно заметить, что типы созревания, в которых встречаются выбросы по содержанию белка и масла совпадают (в групах 2,3,4 и 5). Кроме того, можно заметить, что выбросы, которые детектируются по большей части обратно-пропорциональны друг другу. Логично, что при высоком содержании масла будет низкое содержание белка. Поэтому мы решили не выбрасывать эти значения.
outliers_p <- boxplot.stats(data_without_na$protein_content)$out
plot1 <- ggplot(data_without_na, aes(x = maturation_group, y = protein_content, group = maturation_group)) +
geom_boxplot(
outlier.shape = NA,
color = "black",
lwd = 1.0
) +
geom_jitter(
aes(color = ifelse(protein_content %in% outliers_p, "outlier", "normal")),
width = 0.2,
height = 0,
size = 1
) +
scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +
theme_minimal(base_size = 12) +
labs(
x = "Группа созревания",
y = "Содержание белка %",
) +
theme(
legend.position = "none",
)
outliers_o <- boxplot.stats(data_without_na$oil_content)$out
plot2 <- ggplot(data_without_na, aes(x = maturation_group, y = oil_content, group = maturation_group)) +
geom_boxplot(
outlier.shape = NA,
color = "black",
lwd = 1.0
) +
geom_jitter(
aes(color = ifelse(oil_content %in% outliers_o, "outlier", "normal")),
width = 0.2,
height = 0,
size = 1
) +
scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +
theme_minimal(base_size = 12) +
labs(
x = "Группа созревания",
y = "Содержание масла %",
) +
theme(
legend.position = "none",
)
combined_plot <- (plot1 | plot2) +
plot_annotation(
title = "Поиск выбросов",
theme = theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
)
)
combined_plot
3. Более изощренный поиск выбросов с помощью расстояния махаланобиса.
Давайте попробуем проанализировать все численные данные по годам.
Предположим, что порог для выбросов - 90%. Давайте посмотрим как
распределяются все наши численные данные внутри каждого года. Будем
смотреть в зависимости от года, поскольку в зависимости от погодных
условий в течение года могут меняться и физиологические параметры
растений.
# Визуализация многомерных выбросов по годам
library(ggplot2)
ggplot(numeric_data, aes(x = year, y = mahalanobis_dist, color = outlier)) +
geom_boxplot(outlier.shape = NA, fill = "grey", alpha = 0.5) +
geom_jitter(width = 0.2, height = 0, aes(color = outlier), size = 2) +
scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +
theme_minimal(base_size = 10) +
labs(
title = "Многомерные выбросы по годам",
x = "Год",
y = "Mahalanobis Distance",
color = "Тип выброса"
) +
theme(legend.position = "right")
И несмотря на такой изощренный способ поиска, мы настаиваем, что лучше
оставить значения, которые сильно отклоняются, потому что эти данные
могут говорить о каких-то довольно ценных сортах.
Проверим распределения численных величин. С помощью функции
map в пакете purr возможно значительно
сократить описание функции и переменных, применив
map(var, func)
## Warning: Removed 1543 rows containing non-finite outside the scale range
## (`stat_density()`).
## Warning: Removed 1686 rows containing non-finite outside the scale range
## (`stat_density()`).
## Removed 1686 rows containing non-finite outside the scale range
## (`stat_density()`).
В целом можно отметить, что численные величины распределены относительно равномерно (имееются несклько пиков), но не наблюдаем значительных отклонений.
Знания биологии позволяют нам сказать, что содержание масла прямопропорционально связано с продуктивностью, тогда как содержание белка - обратно. Давайте посмотрим, насколько это воспроизводится на наших данных.
plot_1 <- ggplot(data_without_na, aes(x=productivity, y=oil_content))+
geom_point(shape = 21, size = 1, fill = "#c45824", color = "black") +
theme_minimal() +
labs(
x = "Продуктивность (г/м²)",
y = "Содержание масла (%)"
)
plot_2 <- ggplot(data_without_na, aes(x=productivity, y=protein_content))+
geom_point(shape = 21, size = 1, fill = "#c45824", color = "black") +
theme_minimal() +
labs(
x = "Продуктивность (г/м²)",
y = "Содержание белка (%)"
)
combined_plot <- (plot_1 | plot_2) + plot_annotation(
title = "Распределение данных по содержанию масла и белка",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
)
)
combined_plotПроверим статистичекую значимость кореляций продуктивности с различными параметрами:
Сначала посмотрим на нормальность наших выборок с помощью
shapiro.test
Нулевая гиполеза H0 - выборки распределены нормально
Альтернативная H1 - выборки распределены ненормально
shapiro_var1 <- shapiro.test(data$productivity)
print(paste("p-value для productivity:", shapiro_var1$p.value))## [1] "p-value для productivity: 0.00644333565970214"
shapiro_var2 <- shapiro.test(data$oil_content)
print(paste("p-value для oil_content:", shapiro_var2$p.value))## [1] "p-value для oil_content: 1.36929282750539e-11"
shapiro_var2 <- shapiro.test(data$protein_content)
print(paste("p-value для protein_content:", shapiro_var2$p.value))## [1] "p-value для protein_content: 1.38191116472042e-06"
Видим, что выборки распределены ненормально (p-value<0.5), применим корреляцию по Спирману (продуктивность и содержание масла):
cor_test_result <- cor.test(data$productivity, data$oil_content, method = "spearman")
print(paste("Коэффициент корреляции Спирмена:", cor_test_result$estimate))## [1] "Коэффициент корреляции Спирмена: 0.514123548211428"
## [1] "p-value: 1.53882917165075e-85"
Продуктивность и содержание белков:
cor_test_result <- cor.test(data$productivity, data$protein_content, method = "spearman")
print(paste("Коэффициент корреляции Спирмена:", cor_test_result$estimate))## [1] "Коэффициент корреляции Спирмена: -0.382382602249352"
## [1] "p-value: 6.19818358642783e-45"
Согласно шкале Чедокка, при абсолютном значении коэффициента корреляции Спирмена < 0.03 корреляция слабая, тогда как при от 0.5 до 0.7 заметная. В нашем случае мы обнаружили отрицательную корреляцию между продуктивностью и содержанием белков, тогда как продуктивность и содержание масел коррелируют более значимо.
Взглянем на коэфициенты, используя для этого пакет
psych.
Прежде отберем в отдельный датафрейм интересующие нас количественные переменные
## productivity oil_content protein_content
## productivity 1.0000000 0.5475153 -0.4543480
## oil_content 0.5475153 1.0000000 -0.5929773
## protein_content -0.4543480 -0.5929773 1.0000000
Наблюдаем отицательную корреляцию между продуктивностью, группой цветения и содержанием белков. Корреляция положительна между продуктивностью и содержанием масла. Всё как по учебнику.
H0 - место происхождения не влияет на продуктивность
H1 - место происхождения влияет на продуктивность
Посмотрим на боксплоты продуктивности в зависимости от места происхождения:
По построенным боксплотам видно, что среднее групп различается, поэтому
гипотезу стоит проверить
## Df Sum Sq Mean Sq F value Pr(>F)
## origin 27 792464 29351 8.771 <2e-16 ***
## Residuals 1167 3905291 3346
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
По полученным графикам построенной модели, определим соблюдение допущений:
Остатки распределены нормально и линейно - допущение выполняется, однако есть проблемы с графиком дисперсий, попробуем численные тесты для проверки:
## Warning: Variances differ between groups (Bartlett Test, p = 0.000).
p-value < 0.05, значит, что дисперсии между группами не равны.
Попробуем провести тест анова Уэлча, которая является модификацией теста классической анова не требующей равенства дисперсий:
##
## One-way analysis of means (not assuming equal variances)
##
## data: productivity and origin
## F = 85.125, num df = 27.000, denom df = 31.038, p-value < 2.2e-16
p-value - значительно ниже 0.05, поэтому мы можем отклонить нулевую гипотезу и предположить, что место происхождения влияет на продуктивность сортов.
Однако теперь невозможно использовать классические post hoc тесты, потому что они необходимы для моделей с равенством дисперсий. Все-таки попробуем создать модель несмотря на разницу дисперсий между группами и уже потом провести poc-hoc тесты.
Теперь попробуем выяснить, между каким группами возникают отличия. Для этого проведем тест Тьюки (Tukey HSD), который более строгий и выделяет следующие отличия между странами:
Тест Тьюки более строгий и выделяет следующие отличия между странами:
China - Austria: p = 0.0009
Japan - Austria: p = 0.0010
USA - Austria: p = 0.0092
China - Belarus: p = 0.0013
Japan - Belarus: p = 0.0009
Russia - China: p = 0.0000
Ukraine - China: p = 0.0035
Czech - China: p = 0.0110
Japan - France: p = 0.0216
USA - Czech: p = 0.0223
Отметим, что отличия между местами произрастания замечены в странах находящихся в разных частях света и соотвественно различных климатических зонах: можно выделить страны СНГ (Россия, Беларусь, Украина), страны Азии (Китай, Япония), страны европейской части (Чехия, Франция, Австрия) и Америку.
При поиске выбросов мы предполагали, что гетерогенность выборки достигается не только сортовыми различиями, но и тем, что растения в разные года различались климатические условия, что несомненно влияло на растения. Давайте проверим это гипотезу.
H0 - продуктивность не менялась от года к году
H1 - продуктивность отличалась в разные года
Сначала визуально оценим распределение
По построенным боксплотам видно, что среднее групп различается, поэтому гипотезу стоит проверить
## Df Sum Sq Mean Sq F value Pr(>F)
## year 4 351522 87880 24.06 <2e-16 ***
## Residuals 1190 4346234 3652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
По полученным графикам построенной модели, определим соблюдение допущений:
Полученные графики кажутся приемлемыми для использования модели и проведения post hoc теста.
В данном случае решили проверить в каких граппах есть значимые отличия с помощью SNK теста:
plot(agricolae::SNK.test(model1, trt = c('year')), main = "Результаты теста SNK по годам")
mtext("Продуктивность", side = 2, line = 3) # Подпись оси Y
mtext("Год", side = 1, line = 2)Наиболее сильно климатические условия на урожайность повлияли 2019 и в 2021, также образовалась отдельная группа, внутри которой статистически значимых влияний на продуктивность не обнаружено, это: 2018, 2017, 2020 года.
Возьмём все числовые переменные в качестве предикторов для предсказания продуктивности
model_prod <- lm(productivity ~ flowering_group+maturation_group+vegetation_period+protein_content+oil_content, data_without_na)
summary(model_prod)##
## Call:
## lm(formula = productivity ~ flowering_group + maturation_group +
## vegetation_period + protein_content + oil_content, data = data_without_na)
##
## Residuals:
## Min 1Q Median 3Q Max
## -137.041 -31.691 -0.026 31.938 177.910
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -135.2531 46.7716 -2.892 0.003901 **
## flowering_group.L -49.7753 13.5435 -3.675 0.000248 ***
## flowering_group.Q -1.8186 10.1731 -0.179 0.858152
## flowering_group.C -6.8416 8.1239 -0.842 0.399873
## flowering_group^4 -11.5272 8.1703 -1.411 0.158546
## flowering_group^5 5.5259 10.4461 0.529 0.596913
## flowering_group^6 -6.3540 6.3646 -0.998 0.318327
## flowering_group^7 0.0156 9.9404 0.002 0.998748
## flowering_group^8 8.6407 9.3895 0.920 0.357631
## maturation_group.L -33.0893 9.3378 -3.544 0.000410 ***
## maturation_group.Q -45.1036 6.7629 -6.669 3.94e-11 ***
## maturation_group.C 6.4134 4.4651 1.436 0.151172
## maturation_group^4 -6.7338 2.7922 -2.412 0.016031 *
## vegetation_period 1.9783 0.1654 11.958 < 2e-16 ***
## protein_content -4.2524 0.6458 -6.584 6.86e-11 ***
## oil_content 12.6826 1.0111 12.544 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 45.86 on 1179 degrees of freedom
## Multiple R-squared: 0.4721, Adjusted R-squared: 0.4653
## F-statistic: 70.28 on 15 and 1179 DF, p-value: < 2.2e-16
Согласно значениям Pr(>|t|) все используемые предикторы являются важными. Но в любом случае попробуем проанализировать как измениться предсказательность модели при исключении одного/нескольких предикторов при обнаружении мультиколлинеарности между ними.
Посмотрим на графиках
Оценим нормальность распределения остатков тестом Шапиро-Уилка
##
## Shapiro-Wilk normality test
##
## data: residuals(model_prod)
## W = 0.99869, p-value = 0.53
p-value = 0.1262 - остатки подчиняются нормальному распределению
Оценим коллинеарность предикторов с помощью
corr.test
## vegetation_period protein_content oil_content
## vegetation_period 1 NA NA
## protein_content NA 1 NA
## oil_content NA NA 1
Многие из предикторов сильно коррелируют, например группа цветения c вегетационным периодом и группой созревания.
Оценим с помощью коэффициента инфляции дисперсии (VIF) из пакета car мультиколлениарность предикторов:
## GVIF Df GVIF^(1/(2*Df))
## flowering_group 8.397932 8 1.142249
## maturation_group 10.424573 4 1.340471
## vegetation_period 2.223400 1 1.491107
## protein_content 1.618014 1 1.272012
## oil_content 1.718886 1 1.311063
Получили значения, не превышающие 5, однако попробуем убрать из ряда предикторов с наибольшим значением. Например, такие, как maturation_group и vegetation_period. Начнем с maturation_group
## GVIF Df GVIF^(1/(2*Df))
## flowering_group 1.912530 8 1.041359
## vegetation_period 1.669979 1 1.292277
## protein_content 1.560378 1 1.249151
## oil_content 1.677473 1 1.295173
Значения коэффициентов vif немного снизились для остальных предикторов.
Попробуем оценить значимость всех предикторов с помощью drop1
## Single term deletions
##
## Model:
## productivity ~ flowering_group + maturation_group + vegetation_period +
## protein_content + oil_content
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 2480123 9159.3
## flowering_group 8 64676 2544799 9174.1 3.8432 0.0001775 ***
## maturation_group 4 202289 2682412 9245.0 24.0410 < 2.2e-16 ***
## vegetation_period 1 300804 2780926 9294.1 142.9959 < 2.2e-16 ***
## protein_content 1 91197 2571320 9200.5 43.3531 6.865e-11 ***
## oil_content 1 330983 2811106 9307.0 157.3424 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Согласно оценке, наибольшую значимость имеют предикторы: oil_content, vegetation_period. Наименее значимый - группа цветения (flowering_group).
Сравним модели с помощью информационного критерия Акаике (AIC), Adjusted R-squared и выберем наиболее подходящую.
##
## Call:
## lm(formula = productivity ~ flowering_group + vegetation_period +
## protein_content + oil_content, data = data_without_na)
##
## Residuals:
## Min 1Q Median 3Q Max
## -137.558 -32.904 1.545 32.888 162.343
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -76.9648 46.3751 -1.660 0.0973 .
## flowering_group.L -57.5438 9.8652 -5.833 7.02e-09 ***
## flowering_group.Q -37.5436 8.9469 -4.196 2.92e-05 ***
## flowering_group.C -4.9303 8.3161 -0.593 0.5534
## flowering_group^4 -1.0774 8.1919 -0.132 0.8954
## flowering_group^5 -2.7247 10.5309 -0.259 0.7959
## flowering_group^6 3.9505 6.0728 0.651 0.5155
## flowering_group^7 -0.2820 9.8257 -0.029 0.9771
## flowering_group^8 16.5387 9.4797 1.745 0.0813 .
## vegetation_period 1.4209 0.1489 9.545 < 2e-16 ***
## protein_content -4.6166 0.6585 -7.011 3.97e-12 ***
## oil_content 13.3814 1.0370 12.904 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 47.62 on 1183 degrees of freedom
## Multiple R-squared: 0.429, Adjusted R-squared: 0.4237
## F-statistic: 80.8 on 11 and 1183 DF, p-value: < 2.2e-16
## Single term deletions
##
## Model:
## productivity ~ flowering_group + vegetation_period + protein_content +
## oil_content
## Df Sum of Sq RSS AIC
## <none> 2682412 9245.0
## flowering_group 8 367004 3049416 9382.2
## vegetation_period 1 206602 2889014 9331.7
## protein_content 1 111457 2793869 9291.7
## oil_content 1 377556 3059968 9400.4
Для модели 1 имеем:
Adjusted R-squared: 0.3918
AIC: 9302.5
Модель с исключением из предикторов группы цветения (flowering_group)
##
## Call:
## lm(formula = productivity ~ maturation_group + vegetation_period +
## protein_content + oil_content, data = data_without_na)
##
## Residuals:
## Min 1Q Median 3Q Max
## -136.789 -32.833 0.273 31.376 166.425
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -133.0844 46.4649 -2.864 0.00425 **
## maturation_group.L -57.3939 6.8717 -8.352 < 2e-16 ***
## maturation_group.Q -46.9199 4.5602 -10.289 < 2e-16 ***
## maturation_group.C 0.3291 3.5852 0.092 0.92687
## maturation_group^4 -4.7701 2.6328 -1.812 0.07027 .
## vegetation_period 1.8183 0.1628 11.169 < 2e-16 ***
## protein_content -4.1793 0.6484 -6.446 1.67e-10 ***
## oil_content 13.2361 0.9973 13.271 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 46.3 on 1187 degrees of freedom
## Multiple R-squared: 0.4583, Adjusted R-squared: 0.4551
## F-statistic: 143.5 on 7 and 1187 DF, p-value: < 2.2e-16
## Single term deletions
##
## Model:
## productivity ~ maturation_group + vegetation_period + protein_content +
## oil_content
## Df Sum of Sq RSS AIC
## <none> 2544799 9174.1
## maturation_group 4 504616 3049416 9382.2
## vegetation_period 1 267435 2812234 9291.5
## protein_content 1 89077 2633876 9213.2
## oil_content 1 377609 2922408 9337.4Для модели 2 имеем:
Adjusted R-squared: 0.4033
AIC: 9279.7
Модель с исключением из предикторов данных о вегетационном периоде (vegetation_period)
##
## Call:
## lm(formula = productivity ~ flowering_group + maturation_group +
## protein_content + oil_content, data = data_without_na)
##
## Residuals:
## Min 1Q Median 3Q Max
## -153.765 -32.527 1.365 34.225 164.032
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.4945 44.8931 2.239 0.0254 *
## flowering_group.L -33.0707 14.2588 -2.319 0.0205 *
## flowering_group.Q -2.5405 10.7676 -0.236 0.8135
## flowering_group.C -5.7597 8.5983 -0.670 0.5031
## flowering_group^4 -12.7773 8.6472 -1.478 0.1398
## flowering_group^5 7.3222 11.0556 0.662 0.5079
## flowering_group^6 -8.4767 6.7341 -1.259 0.2084
## flowering_group^7 7.3414 10.5015 0.699 0.4846
## flowering_group^8 2.5783 9.9239 0.260 0.7951
## maturation_group.L 13.7320 8.9728 1.530 0.1262
## maturation_group.Q -46.7610 7.1567 -6.534 9.51e-11 ***
## maturation_group.C 9.6968 4.7172 2.056 0.0400 *
## maturation_group^4 -6.3021 2.9551 -2.133 0.0332 *
## protein_content -5.1364 0.6791 -7.563 7.87e-14 ***
## oil_content 12.9425 1.0699 12.096 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 48.55 on 1180 degrees of freedom
## Multiple R-squared: 0.408, Adjusted R-squared: 0.401
## F-statistic: 58.1 on 14 and 1180 DF, p-value: < 2.2e-16
## Single term deletions
##
## Model:
## productivity ~ flowering_group + maturation_group + protein_content +
## oil_content
## Df Sum of Sq RSS AIC
## <none> 2780926 9294.1
## flowering_group 8 31308 2812234 9291.5
## maturation_group 4 108088 2889014 9331.7
## protein_content 1 134817 2915743 9348.7
## oil_content 1 344846 3125772 9431.8Для модели 3 имеем:
Adjusted R-squared: 0.3512
AIC: 9379.7
При этом для исходной модели с сохранением всех предиктров:
Adjusted R-squared: 0.4141
AIC: 9258.8
Таким образом, точность предсказания исходной модели максимальная в сравнении с тремя рассмотренными выше и составляет 41.41%; также для нее информационный критерий Акаике был минимален
Для предсказания продуктивности сортов наиболее подходящей будет модель, учитывающая группу цветения, вегетационный период, группу созревания, содержание белков и масел
model5 <- lme4::lmer(productivity ~ leaf_shape + maturation_group + lodging_type + (1 | year), data = data_without_na)
summary(model5)## Linear mixed model fit by REML ['lmerMod']
## Formula: productivity ~ leaf_shape + maturation_group + lodging_type +
## (1 | year)
## Data: data_without_na
##
## REML criterion at convergence: 12927.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.00933 -0.69312 -0.01955 0.68005 2.88914
##
## Random effects:
## Groups Name Variance Std.Dev.
## year (Intercept) 476.5 21.83
## Residual 3004.9 54.82
## Number of obs: 1195, groups: year, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 151.888 12.917 11.759
## leaf_shaperound 6.679 5.703 1.171
## maturation_group.L -1.175 6.025 -0.195
## maturation_group.Q -71.319 5.231 -13.635
## maturation_group.C 12.295 4.244 2.897
## maturation_group^4 -6.904 3.114 -2.217
## lodging_typeno -4.654 6.832 -0.681
## lodging_typeyes -50.538 12.476 -4.051
##
## Correlation of Fixed Effects:
## (Intr) lf_shp mtr_.L mtr_.Q mtr_.C mtr_^4 ldgng_typn
## leaf_shprnd -0.381
## mtrtn_grp.L -0.059 0.001
## mtrtn_grp.Q 0.155 -0.069 -0.299
## mtrtn_grp.C -0.105 0.129 0.477 -0.221
## mtrtn_grp^4 0.077 -0.092 -0.116 0.381 -0.090
## lodgng_typn -0.481 -0.054 0.035 -0.041 0.047 -0.042
## ldgng_typys -0.279 0.002 -0.097 -0.103 -0.045 -0.057 0.512
## [[1]]
## `geom_smooth()` using formula = 'y ~ x'
##
## [[2]]
## [[2]]$year
## `geom_smooth()` using formula = 'y ~ x'
##
##
## [[3]]
##
## [[4]]
## `geom_smooth()` using formula = 'y ~ x'
Scaled residuals показывают распределение остатков
модели. Нормальное распределение остатков позволяет предположить, что
модель адекватна.
Случайные эффекты:
year имеет вариацию 464.3 и стандартное отклонение
21.55, что говорит о значительном влиянии года на продуктивность.
Фиксированные эффекты:
Intercept: Оценка пересечения равна 202.435. Это
базовый уровень продуктивности для группы, когда все остальные
переменные равны нулю.
leaf_shaperound: Небольшой положительный эффект
(+2.670), но он не статистически значим (t value = 0.439).
maturation_group: Отрицательный значимый эффект
(-8.021), что указывает на то, что увеличение группы созреваемости
связано с уменьшением продуктивности.
lodging_typeno: Отрицательный эффект (-7.796) при
отсутствии уклона, но не статистически значимый.
lodging_typeyes: Заметный отрицательный эффект
(-67.445) с высоким t-значением (-5.053), что указывает на значительное
снижение продуктивности, когда растения предрасположены к
уклону.
Выводы из модели показывают, что maturation_group и
lodging_type имеют значительное влияние на продуктивность.
Возможная предрасположенность к уклону существенно снижает
продуктивность.
1. С помощью дисперсионного анализа ANOVA мы установили, что регион происхождения влияет на продуктивность: отличия между местами произрастания замечены в странах находящихся в разных частях света и соотвественно различных климатических зонах.
2. С помощью дисперсионного анализа ANOVA мы установили, что продуктивность сои также отличается в зависимости от года, с чем связана часть гетерогенности и выбросов в нашей выборке. Наиболее сильно климатические условия на урожайность повлияли 2019 и в 2021 году.
3. Подобрали модель предсказывающую продуктивность которая в качестве предикторов учитывает группу цветения, вегетационный период, группу созревания, содержание белков и масел.